DIPGEOPR_02

class: inverse,left, middle
background-image: url(data:image/png;base64,#background.png)
background-size: cover

##Módulo 2: Análisis espacial y Percepción Remota satelital óptica

###Introducción al análisis masivo: Paralelización ráster

José A. Lastra <a href="http://github.com/JoseLastra">
Github: JoseLastra</a> 
<a href="mailto:jose.lastra@pucv.cl"> jose.lastra@pucv.cl</a> | <a href="mailto:jose.lastramunoz@wur.nl"> jose.lastramunoz@wur.nl</a>

.large[<a href="https://www.pucv.cl/uuaa/site/edic/base/port/labgrs.html">LabGRS</a> | Octubre, 2024]

---
class: center,middle
background-image: url(data:image/png;base64,#labgrs_logo.png)
background-size: 35%

---
## Contenidos

.pull-left[
- Paralelización:
 * Librería parallel
 * Librería foreach y doParallel
 
- Librería terra:
 * app()
 * sapp()
 * lapply()
 ]
 
.pull-right[
<img src="data:image/png;base64,#https://raw.githubusercontent.com/allisonhorst/stats-illustrations/main/rstats-artwork/r_rollercoaster.png" width="650px"/>
© Allison Horst
]

---
## Recordando algunos conceptos

--
- Problemas de cálculo ([Jones, 2017](https://nceas.github.io/oss-lessons/parallel-computing-in-r/parallel-computing-in-r.html)):
  * **CPU**: proceso toma mucho tiempo de CPU.
  * **Memoria**: Demanda demasiada memoria.
  * **I/O**: toma mucho tiempo la lectura/escritura desde el disco.
  * **Network**: La red toma demasiado tiempo en transferir.

--
- La paralelización viene a "solucionar" los problemas asociados a la CPU, dada la presencia de múltiples núcleos en los pc's modernos.

---
## Paralelización: Conceptos

<center><img src="data:image/png;base64,#sock_fork.png" width="900px"/></center>
<center>Fuente: Benito, 2021</center>

---
## foreach y doParallel: Ejemplo práctico con rásters

--
- Crearemos una lista de archivos **MOD13Q1** correspondientes al índice NDVI desde el año 2000 hasta el 2021.

--
- Ejecutaremos dos tareas simples: (1) cortar y (2) enmascarar datos a un área de estudio.

--
- Usaremos los siguientes archivos: (1) **san_javier.gpkg**, (2) **MOD13Q1_dates_ts492.csv** y (3) **MOD13Q1_bands (carpeta)**

``` r
## Listando archivos
bandas <- list.files(path = "MOD13Q1_bands/", pattern = "*.tif", full.names = T)

## archivo de fechas
tabla <- read_csv("MOD13Q1_dates_ts492.csv")

## shape de San Javier
shp <- vect("san_javier.gpkg")
```

---
## foreach y doParallel: Ejemplo práctico con rásters

--
- Preparación de clúster de procesamiento.

``` r
## nombres de salida archivos
nombres <- paste0(
 "San_Javier_MOD13Q1/",
 tabla$dates, "_",
 tabla$ID, "_San_Javier.tif"
)

## Cluster set up
nCluster <- detectCores(logical = F) - 1 # Número de núcleos físicos
cl <- makeCluster(nCluster, type = "PSOCK") # Configuración del cluster
registerDoParallel(cl) # registrar para doParallel
```

---
## foreach y doParallel: Ejemplo práctico con raster

--
- Configuración de foreach y detención del clúster

``` r
foreach(i = 1:length(nombres)) %dopar% {
  # importar librerías a los nodos
  library(raster)
  library(magrittr)

r <- bandas[i] %>% raster() # crear banda en R
 m <- r %>%
 crop(shp) %>%
 mask(shp) # #cortar y enmascarar

# guardar resultado
  writeRaster(m,
    filename = nombres[i],
    overwrite = T, datatype = "INT2S"
  )
}

stopCluster(cl)
```

--
- ¿Qué falló?

---
## foreach y doParallel: Ejemplo práctico con raster

--
- Esto lo podemos solucionar de varias maneras:
 - Utilizar __sf__ para leer nuestro vector: `shp <- read_sf('archivo.gpkg')`
 - Pasar la ruta del archivo y leerlo en cada procesamiento: `ruta <- 'san_javier.gpkg'` y `shp <- vect(ruta)` dentro del __foreach__.
 - Empaquetar nuestro SpatVector `shp <- vect('archivo.gpkg') %>% wrap()` y `vect(shp)` ó `unwrap(shp)` en nuestro __foreach__

---

<center><img src="data:image/png;base64,#cpu.png" width="600px"/></center>
<center>LabGRS, 2024</center>

---
## terra y la paralelización

--
- Los ejemplos de paralelización vistos anteriormente, donde se envía información a cada nodo, no son compatibles con los objetos de la librería **terra** ([ver más Hijmans et al., 2024 pág. 8](https://cran.r-project.org/web/packages/terra/terra.pdf))

--
- terra provee funciones que permiten realizar operaciones empleando múltiples núcleos. Siendo **app()** la más común y simple de emplear.

--
- Podemos ajustar este comportamiento; pero algunas cosas no pueden ser trabajadas de esta forma.

--
- Por esta razón, terra provee funciones que permiten realizar operaciones empleando múltiples núcleos. Siendo __app()__ la más común y directa.

``` r
ndvi_sj <- list.files(
 path = "San_Javier_MOD13Q1/",
 pattern = "*.tif",
 full.names = T
) %>% rast()
```

---
## app()

--
- Esta función nos permite operar en un **SpatRaster** multi o monobanda para aplicar una función normalmente de resumen (sum, min, max, median, etc.). También podemos operar en **SpatRasterDataset**.

--
- **Importante**: dado la eficiencia computacional de **terra** muchas funciones resumen pueden ser aplicadas directamente con un funcionamiento óptimo.

--
- Sin embargo, para archivos muy grandes puede ser útil la aplicación de **app()**

``` r
mean_noApp <- ndvi_sj %>% mean(na.rm = T) # simple

mean_app <- app(ndvi_sj, fun = mean, na.rm = T, cores = 5) # 5 cores
```

--
- **app()** opera los archivos de la misma forma que **apply()**: donde cada banda es una columna en la matriz y los píxeles son las filas.

---

---
## sapp()

--
- Si lo que queremos es aplicar una función a cada banda dentro de nuestro **SpatRaster**, podemos emplear *sapp()*.

``` r
ndvi <- list.files(
 path = "MOD13Q1_bands/",
 pattern = "*.tif",
 full.names = T
) %>% rast()

shp_v <- vect('san_javier.gpkg')

# versión simple
ndvi_sj <- ndvi %>%
 crop(shp_v) %>%
 mask(shp_v)
```

---
## sapp()

``` r
# versión con sapp()
## Creando la función
crop_mask <- function(x, ...) {
 r <- x %>%
 crop(shp_v) %>%
 mask(shp_v)
 r
}

## aplicación
ndvi_sj_sapp <- sapp(ndvi_sj, fun = crop_mask)
```

---

---
## Observaciones

--
- **sapp()** invoca a **sapply()** por detrás y aún no se ha implementado la paralelización para esta alternativa (**Información actualizada al 14 de Octubre de 2024**).

--
- Tal cuál lo indica la documentación: *"(...) La paralelización es soportada, pero a menudo no ayuda, e incluso puede ser más lenta."*

--
- Justificaciones para paralelizar con terra:
  * Disponemos de más de 8 núcleos para procesar.
  * Tenemos un archivo muy grande que trabajar.
  * O tenemos una función muy compleja (lenta) para aplicar sobre nuestros archivos.
  * Muchas veces el uso de **lapply()** puede ser suficiente.

---
## lapply() con terra

--
- Tomemos nuestra lista de archivos originales y la función creada anteriormente para cortar y enmascarar.

``` r
ndvi_list <- list.files(
 path = "MOD13Q1_bands/",
 pattern = "*.tif",
 full.names = T
)

# adaptando la función para lapply
crop_mask <- function(x, ...) {
 x <- rast(x) # SpatRaster
 r <- x %>%
 crop(shp_v) %>%
 mask(shp_v) # aplicación de operación
 rm(x) # liberar memoria
 r # entregar objeto
}

ndvi_crop <- lapply(ndvi_list, FUN = crop_mask)
```

---
## Recomendaciones

--
- Siempre es mejor probar primero el código en un área pequeña antes de ejecutar el procesamiento de gran escala.

--
- Podemos usar un área de prueba transformada en un tibble e ir porbando pixel a pixel la función, facilitando la búsqueda de errores.

--
- También podemos emplear otras estructuras de código como transformar nuestros objetos en listas y data frames para paralelizar empleando __future_lapply()__

---
## Funciones complejas en app()

--
- El principal pro de __terra__ y `app()` tiene que ver con la potencia y eficiencia computacional.

--
- Especialmente, si lo comparamos con su antecesor de librería __raster__: `calc()` y `clusterR()`

--
- Calculemos la tendencia de nuestro NDVI para el SpatRaster __ndvi_crop__

``` r
## Creación de función para evaluar tendencia
## definiendo tiempo
time <- 1:nlyr(ndvi_crop)

calc_trend <- function(x, time) { 
 if(all(is.na(x))){ # Si la serie es solo NA's devolvemos NA
 return(rep(NA,2))
 } else { 
 m <- lm(x ~ time)$coefficients[2] # extraccion pendiente
 p <- summary(lm(x ~ time))$coefficients[2,4] # extraccion pvalue
 out_vals <- c(m,p) # creacion de salida
 names(out_vals) <- c('trend', 'pval') # nombres para identificar variables
 return(out_vals)
 }
}
```

---
---
## Funciones complejas en app()

``` r
## aplicacion de funcion al archivo spatRaster
NDVI_trend <- ndvi_crop %>% 
 app(fun = calc_trend, time = time, cores = 6) 
```

---
## Bibliografía complementaria

- Bosco, J. (2018). "R Para Principantes". [Ver](https://bookdown.org/jboscomendoza/r-principiantes4/)

- Cao, R. y Fernández, R. (2022). "Introducción al procesamiento en paralelo en R". En: "Técnicas de Remuestreo". [Ver](https://rubenfcasal.github.io/book_remuestreo/intro-hpc.html)

- Gillespie, C., & Lovelace, R. (2021). "Efficient R programming: a practical guide to smarter programming." O'Reilly Media, Inc. [Ver](https://csgillespie.github.io/efficientR/)

- Hijmans, R. J., Bivand, R., Forner, K., Ooms, J., Pebesma, E., & Sumner, M. D. (2022). Package ‘terra’. [Ver](https://brieger.esalq.usp.br/CRAN/web/packages/terra/terra.pdf)

- Jones, M. (2017). "Quick Intro to Parallel Computing in R". [Ver](https://nceas.github.io/oss-lessons/parallel-computing-in-r/parallel-computing-in-r.html)

- McCallum, E., & Weston, S. (2011). "Parallel R". O'Reilly Media, Inc.

- Orellana, J. (2018). "HPC con R para investigadores". [Ver](https://bookdown.org/content/1498/)

---
class: middle
![logo](data:image/png;base64,#LABGRS-PUCV.png)